R Notebook

The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file

Setup

Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here

#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ─────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))

processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd

  |                                                                                                         
  |                                                                                                   |   0%
  |                                                                                                         
  |....                                                                                               |   4%
  |                                                                                                         
  |.........                                                                                          |   9%
  |                                                                                                         
  |.............                                                                                      |  13%
  |                                                                                                         
  |.................                                                                                  |  17%
  |                                                                                                         
  |......................                                                                             |  22%
  |                                                                                                         
  |..........................                                                                         |  26%
  |                                                                                                         
  |..............................                                                                     |  30%
  |                                                                                                         
  |..................................                                                                 |  35%
  |                                                                                                         
  |.......................................                                                            |  39%
  |                                                                                                         
  |...........................................                                                        |  43%
  |                                                                                                         
  |...............................................                                                    |  48%
  |                                                                                                         
  |....................................................                                               |  52%
  |                                                                                                         
  |........................................................                                           |  57%
  |                                                                                                         
  |............................................................                                       |  61%
  |                                                                                                         
  |.................................................................                                  |  65%
  |                                                                                                         
  |.....................................................................                              |  70%
  |                                                                                                         
  |.........................................................................                          |  74%
  |                                                                                                         
  |.............................................................................                      |  78%
  |                                                                                                         
  |..................................................................................                 |  83%
  |                                                                                                         
  |......................................................................................             |  87%
  |                                                                                                         
  |..........................................................................................         |  91%
  |                                                                                                         
  |...............................................................................................    |  96%
  |                                                                                                         
  |...................................................................................................| 100%

output file: /tmp/RtmpjHvOEV/file136755260521

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Attaching package: ‘flextable’

The following object is masked from ‘package:purrr’:

    compose


Attaching package: ‘ftExtra’

The following object is masked from ‘package:flextable’:

    separate_header

Warning: Assuming taxa are rows

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

raDf <- nonSpikes_Remake %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes_Remake %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

seqDep <- countMat %>% apply(1, sum)
min(seqDep)

[1] 852

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2 
  4.310891   5.165546   4.689957   5.913745   5.111418   3.870738   5.557505   4.510650   4.848269 
 3-1-S-180   3-1-S-20    3-1-S-5  3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20    3-2-B-5  3-2-B-500 
  4.696624   5.233180   4.817138   4.421646   4.789616   4.646308   5.183631   5.534260   4.990589 
  3-2-B-53  3-2-S-0-2  3-2-S-1-2  3-2-S-180   3-2-S-20    3-2-S-5  3-2-S-500   3-2-S-53  3-3-B-0-2 
  4.976637   3.723508   4.915513   4.789992   4.917679   5.140993   4.974091   4.328113   4.365590 
 3-3-B-1-2  3-3-B-180   3-3-B-20    3-3-B-5  3-3-B-500   3-3-B-53  3-3-S-180   3-3-S-20  3-3-S-500 
  4.883612   3.466347   5.650654   5.357062   5.314725   5.519214   4.983904   4.954011   4.912663 
  3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180   4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53  4-3-O-1-2 
  4.427651   4.278936   4.953905   4.517047   4.315042   4.678940   4.416019   4.274373   5.007986 
 4-3-O-180    4-3-O-5  4-3-O-500   4-3-O-53  4-3-S-0-2  4-3-S-180   4-3-S-20  4-3-S-500   4-3-S-53 
  4.717124   5.223329   4.046068   4.686227   2.805419   4.501172   4.556331   4.744078   4.495019 
 5-1-S-1-2  5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53  5-5-B-0-2  5-5-B-180    5-5-B-5 
  4.457156   4.528074   4.165735   3.886173   4.045319   4.008392   4.584748   5.178673   5.438522 
 5-5-B-500   5-5-B-53  5-5-S-180    5-5-S-5  5-5-S-500   5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 
  5.030883   4.915445   4.358304   4.955525   4.838466   4.253475   4.239569   4.724150   4.847382 
 C_5P1B_20 C_5P1B_500  C_5P1B_53 
  5.376668   4.837291   5.029040 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2 
0.012002661 0.005695199 0.012624379 0.003374655 0.007066477 0.041845812 0.008999821 0.008826596 0.008569315 
  3-1-S-180    3-1-S-20     3-1-S-5   3-2-B-0-2   3-2-B-1-2   3-2-B-180    3-2-B-20     3-2-B-5   3-2-B-500 
0.009531201 0.008122208 0.008588573 0.011376936 0.008519080 0.012306848 0.006157532 0.007217634 0.010603357 
   3-2-B-53   3-2-S-0-2   3-2-S-1-2   3-2-S-180    3-2-S-20     3-2-S-5   3-2-S-500    3-2-S-53   3-3-B-0-2 
0.010199187 0.019991991 0.008528156 0.008806712 0.008659411 0.008230593 0.008635575 0.013009545 0.014582521 
  3-3-B-1-2   3-3-B-180    3-3-B-20     3-3-B-5   3-3-B-500    3-3-B-53   3-3-S-180    3-3-S-20   3-3-S-500 
0.007702219 0.066660513 0.005222085 0.006729824 0.006906726 0.004171370 0.009138070 0.008108037 0.010812581 
   3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180    4-3-B-20     4-3-B-5   4-3-B-500    4-3-B-53   4-3-O-1-2 
0.009967698 0.011272718 0.009273930 0.008735626 0.009611597 0.008821531 0.007439386 0.009714484 0.008645639 
  4-3-O-180     4-3-O-5   4-3-O-500    4-3-O-53   4-3-S-0-2   4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53 
0.010714968 0.008686394 0.015126253 0.010410610 0.124685269 0.014872533 0.010144341 0.009760949 0.012606223 
  5-1-S-1-2   5-1-S-180    5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53   5-5-B-0-2   5-5-B-180     5-5-B-5 
0.011975164 0.012247374 0.013105682 0.018611254 0.016107958 0.014028890 0.008750354 0.008721128 0.006220032 
  5-5-B-500    5-5-B-53   5-5-S-180     5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2  C_5P1B_180  C_5P1B_1P2 
0.008597596 0.009145546 0.013668888 0.010845538 0.009159701 0.014943540 0.008794075 0.007531669 0.005760406 
  C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.004266727 0.007981538 0.005106721 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied observed species numbers

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-222.409  -42.169    3.618   50.716  200.239 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          228869.125 135148.104   1.693 0.094877 .  
log(Size_Class)          32.876      8.188   4.015 0.000149 ***
I(log(Size_Class)^2)     -6.386      1.523  -4.192 8.07e-05 ***
lat                  -11916.986   7042.393  -1.692 0.095123 .  
I(lat^2)                155.218     91.678   1.693 0.094952 .  
depth                     3.348      3.243   1.032 0.305574    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 76.98 on 69 degrees of freedom
Multiple R-squared:  0.2574,    Adjusted R-squared:  0.2035 
F-statistic: 4.782 on 5 and 69 DF,  p-value: 0.0008238

Richness

Rarified chao1 estimates

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-537.58 -134.92  -16.72  110.55 1001.08 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          947118.919 443379.266   2.136 0.036218 *  
log(Size_Class)          85.353     26.863   3.177 0.002224 ** 
I(log(Size_Class)^2)    -17.651      4.998  -3.532 0.000741 ***
lat                  -49369.792  23103.920  -2.137 0.036157 *  
I(lat^2)                643.296    300.769   2.139 0.035990 *  
depth                    19.724     10.640   1.854 0.068034 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 252.5 on 69 degrees of freedom
Multiple R-squared:  0.2128,    Adjusted R-squared:  0.1557 
F-statistic: 3.729 on 5 and 69 DF,  p-value: 0.004785

As above but without latitude and depth

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-468.19 -159.00  -14.57  110.38 1095.72 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           561.748     45.961  12.222  < 2e-16 ***
log(Size_Class)        85.808     27.295   3.144 0.002424 ** 
I(log(Size_Class)^2)  -18.065      5.077  -3.558 0.000666 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 256.9 on 72 degrees of freedom
Multiple R-squared:  0.1498,    Adjusted R-squared:  0.1262 
F-statistic: 6.342 on 2 and 72 DF,  p-value: 0.002904

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3141 -0.1915  0.0272  0.3175  0.7540 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.361e+03  7.921e+02   1.718   0.0903 .  
log(Size_Class)       2.068e-01  4.799e-02   4.310 5.31e-05 ***
I(log(Size_Class)^2) -3.698e-02  8.929e-03  -4.142 9.60e-05 ***
lat                  -7.064e+01  4.127e+01  -1.712   0.0915 .  
I(lat^2)              9.198e-01  5.373e-01   1.712   0.0914 .  
depth                 1.177e-02  1.901e-02   0.619   0.5379    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4512 on 69 degrees of freedom
Multiple R-squared:  0.2981,    Adjusted R-squared:  0.2472 
F-statistic:  5.86 on 5 and 69 DF,  p-value: 0.0001437

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013980 -0.004919 -0.002022  0.000656  0.099590 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.323e+01  2.665e+01  -0.872   0.3865  
log(Size_Class)      -4.000e-03  1.615e-03  -2.478   0.0157 *
I(log(Size_Class)^2)  6.964e-04  3.004e-04   2.318   0.0234 *
lat                   1.209e+00  1.389e+00   0.871   0.3869  
I(lat^2)             -1.572e-02  1.808e-02  -0.870   0.3875  
depth                -3.007e-04  6.395e-04  -0.470   0.6397  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01518 on 69 degrees of freedom
Multiple R-squared:  0.09907,   Adjusted R-squared:  0.03379 
F-statistic: 1.518 on 5 and 69 DF,  p-value: 0.1957

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12 
226.1310 226.1310 200.6366 200.6366 204.4382 161.4112 161.4112 181.9116 200.7873 301.3648 301.3648 275.8704 
      13       14       15       16       17       18       19       20       21       22       23       24 
275.8704 279.6720 236.6450 236.6450 257.1453 257.1453 331.9534 331.9534 306.4590 306.4590 310.2606 267.2336 
      25       26       27       28       29       30       31       32       33       34       35       36 
267.2336 287.7339 306.6096 306.6096 336.7613 336.7613 311.2669 311.2669 315.0685 315.0685 272.0415 272.0415 
      37       38       39       40       41       42       43       44       45       46       47       48 
292.5418 292.5418 325.4489 299.9545 299.9545 303.7561 303.7561 260.7291 260.7291 260.7291 281.2294 281.2294 
      49       50       51       52       53       54       55       56       57       58       59       60 
300.1051 300.1051 294.1015 294.1015 268.6071 268.6071 272.4087 272.4087 229.3817 229.3817 229.3817 249.8820 
      61       62       63       64       65       66       67       68       69       70       71       72 
249.8820 268.7577 268.7577 253.2659 227.7715 227.7715 231.5731 231.5731 188.5461 188.5461 188.5461 209.0464 
      73       74       75 
209.0464 227.9221 227.9221 

$se.fit
 [1] 26.93010 26.93010 26.19587 26.19587 26.05607 27.72062 27.72062 29.27714 33.96388 18.89110 18.89110
[12] 18.27607 18.27607 17.16336 19.99663 19.99663 21.45631 21.45631 19.21214 19.21214 18.75132 18.75132
[23] 17.19863 20.04323 20.04323 21.15902 28.11574 28.11574 19.20997 19.20997 18.71830 18.71830 16.92495
[34] 16.92495 19.56830 19.56830 20.53067 20.53067 18.52166 17.88401 17.88401 15.96344 15.96344 18.40305
[45] 18.40305 18.40305 19.39190 19.39190 26.49688 26.49688 19.27023 19.27023 18.38518 18.38518 16.64424
[56] 16.64424 18.38669 18.38669 18.38669 19.45656 19.45656 26.05823 26.05823 24.40667 23.45563 23.45563
[67] 22.29657 22.29657 23.09072 23.09072 23.09072 24.08798 24.08798 29.16339 29.16339

$df
[1] 69

$residual.scale
[1] 76.97868

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12 
470.9464 470.9464 357.1670 357.1670 439.2520 269.8295 269.8295 392.2784 361.0917 669.0136 669.0136 555.2343 
      13       14       15       16       17       18       19       20       21       22       23       24 
555.2343 637.3192 467.8968 467.8968 590.3456 590.3456 745.6874 745.6874 631.9080 631.9080 713.9930 544.5705 
      25       26       27       28       29       30       31       32       33       34       35       36 
544.5705 667.0194 635.8327 635.8327 751.3244 751.3244 637.5450 637.5450 719.6300 719.6300 550.2076 550.2076 
      37       38       39       40       41       42       43       44       45       46       47       48 
672.6564 672.6564 714.6754 600.8960 600.8960 682.9810 682.9810 513.5586 513.5586 513.5586 636.0074 636.0074 
      49       50       51       52       53       54       55       56       57       58       59       60 
604.8207 604.8207 621.2768 621.2768 507.4974 507.4974 589.5824 589.5824 420.1599 420.1599 420.1599 542.6088 
      61       62       63       64       65       66       67       68       69       70       71       72 
542.6088 511.4221 511.4221 502.7609 388.9815 388.9815 471.0664 471.0664 301.6440 301.6440 301.6440 424.0929 
      73       74       75 
424.0929 392.9061 392.9061 

$se.fit
 [1]  88.34936  88.34936  85.94058  85.94058  85.48195  90.94281  90.94281  96.04928 111.42503  61.97588
[11]  61.97588  59.95815  59.95815  56.30769  65.60278  65.60278  70.39155  70.39155  63.02912  63.02912
[21]  61.51731  61.51731  56.42339  65.75566  65.75566  69.41624  92.23906  92.23906  63.02200  63.02200
[31]  61.40898  61.40898  55.52553  55.52553  64.19755  64.19755  67.35480  67.35480  60.76386  58.67191
[41]  58.67191  52.37114  52.37114  60.37472  60.37472  60.37472  63.61886  63.61886  86.92810  86.92810
[51]  63.21969  63.21969  60.31612  60.31612  54.60461  54.60461  60.32105  60.32105  60.32105  63.83099
[61]  63.83099  85.48902  85.48902  80.07075  76.95068  76.95068  73.14818  73.14818  75.75353  75.75353
[71]  75.75353  79.02524  79.02524  95.67609  95.67609

$df
[1] 69

$residual.scale
[1] 252.5433

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12 
4.483447 4.483447 4.343687 4.343687 4.274961 4.008305 4.008305 4.078398 4.341968 4.948626 4.948626 4.808865 
      13       14       15       16       17       18       19       20       21       22       23       24 
4.808865 4.740140 4.473483 4.473483 4.543576 4.543576 5.149242 5.149242 5.009481 5.009481 4.940756 4.674100 
      25       26       27       28       29       30       31       32       33       34       35       36 
4.674100 4.744193 5.007763 5.007763 5.199876 5.199876 5.060115 5.060115 4.991390 4.991390 4.724733 4.724733 
      37       38       39       40       41       42       43       44       45       46       47       48 
4.794826 4.794826 5.150379 5.010618 5.010618 4.941893 4.941893 4.675236 4.675236 4.675236 4.745329 4.745329 
      49       50       51       52       53       54       55       56       57       58       59       60 
5.008900 5.008900 4.988925 4.988925 4.849165 4.849165 4.780439 4.780439 4.513783 4.513783 4.513783 4.583876 
      61       62       63       64       65       66       67       68       69       70       71       72 
4.583876 4.847446 4.847446 4.769215 4.629455 4.629455 4.560730 4.560730 4.294073 4.294073 4.294073 4.364166 
      73       74       75 
4.364166 4.627736 4.627736 

$se.fit
 [1] 0.15783393 0.15783393 0.15353071 0.15353071 0.15271137 0.16246707 0.16246707 0.17158964 0.19905804
[10] 0.11071836 0.11071836 0.10711375 0.10711375 0.10059229 0.11719773 0.11719773 0.12575275 0.12575275
[19] 0.11259995 0.11259995 0.10989914 0.10989914 0.10079898 0.11747085 0.11747085 0.12401038 0.16478279
[28] 0.16478279 0.11258723 0.11258723 0.10970561 0.10970561 0.09919497 0.09919497 0.11468732 0.11468732
[37] 0.12032767 0.12032767 0.10855312 0.10481590 0.10481590 0.09355973 0.09355973 0.10785794 0.10785794
[46] 0.10785794 0.11365350 0.11365350 0.15529489 0.15529489 0.11294041 0.11294041 0.10775324 0.10775324
[55] 0.09754977 0.09754977 0.10776205 0.10776205 0.10776205 0.11403247 0.11403247 0.15272400 0.15272400
[64] 0.14304440 0.13747047 0.13747047 0.13067740 0.13067740 0.13533179 0.13533179 0.13533179 0.14117663
[73] 0.14117663 0.17092295 0.17092295

$df
[1] 69

$residual.scale
[1] 0.4511623

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9 
0.019539795 0.019539795 0.021890799 0.021890799 0.021592773 0.025095028 0.025095028 0.022774559 0.018846627 
         10          11          12          13          14          15          16          17          18 
0.010590975 0.010590975 0.012941978 0.012941978 0.012643952 0.016146208 0.016146208 0.013825738 0.013825738 
         19          20          21          22          23          24          25          26          27 
0.006662563 0.006662563 0.009013567 0.009013567 0.008715541 0.012217797 0.012217797 0.009897327 0.005969395 
         28          29          30          31          32          33          34          35          36 
0.005969395 0.005562772 0.005562772 0.007913776 0.007913776 0.007615750 0.007615750 0.011118005 0.011118005 
         37          38          39          40          41          42          43          44          45 
0.008797536 0.008797536 0.006391978 0.008742981 0.008742981 0.008444956 0.008444956 0.011947211 0.011947211 
         46          47          48          49          50          51          52          53          54 
0.011947211 0.009626742 0.009626742 0.005698810 0.005698810 0.009303238 0.009303238 0.011654241 0.011654241 
         55          56          57          58          59          60          61          62          63 
0.011356216 0.011356216 0.014858471 0.014858471 0.014858471 0.012538002 0.012538002 0.008610070 0.008610070 
         64          65          66          67          68          69          70          71          72 
0.013332733 0.015683736 0.015683736 0.015385710 0.015385710 0.018887966 0.018887966 0.018887966 0.016567496 
         73          74          75 
0.016567496 0.012639564 0.012639564 

$se.fit
 [1] 0.005310081 0.005310081 0.005165306 0.005165306 0.005137741 0.005465956 0.005465956 0.005772871
 [9] 0.006697004 0.003724950 0.003724950 0.003603678 0.003603678 0.003384274 0.003942939 0.003942939
[17] 0.004230759 0.004230759 0.003788253 0.003788253 0.003697389 0.003697389 0.003391228 0.003952127
[25] 0.003952127 0.004172140 0.005543865 0.005543865 0.003787825 0.003787825 0.003690877 0.003690877
[33] 0.003337263 0.003337263 0.003858479 0.003858479 0.004048241 0.004048241 0.003652104 0.003526371
[41] 0.003526371 0.003147674 0.003147674 0.003628715 0.003628715 0.003628715 0.003823698 0.003823698
[49] 0.005224659 0.005224659 0.003799707 0.003799707 0.003625193 0.003625193 0.003281913 0.003281913
[57] 0.003625490 0.003625490 0.003625490 0.003836448 0.003836448 0.005138166 0.005138166 0.004812510
[65] 0.004624984 0.004624984 0.004396441 0.004396441 0.004553031 0.004553031 0.004553031 0.004749672
[73] 0.004749672 0.005750441 0.005750441

$df
[1] 69

$residual.scale
[1] 0.01517867

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
  bold(i = ~ p< 0.05, j = "p") %>%
  colformat_md() %>%
  set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Metric	Term	Estimate	Standard Error	T Value	p
Observed ASVs	Intercept	2.3 x 105	1.4 x 105	1.69	0.095
	log(Size Class)	3.3 x 101	8.2 x 100	4.01	< 0.001
	log(Size Class)2	-6.4 x 100	1.5 x 100	-4.19	< 0.001
	Latitude	-1.2 x 104	7.0 x 103	-1.69	0.095
	Latitude2	1.6 x 102	9.2 x 101	1.69	0.095
	Depth	3.3 x 100	3.2 x 100	1.03	0.306
Richness (Chao1)	Intercept	9.5 x 105	4.4 x 105	2.14	0.036
	log(Size Class)	8.5 x 101	2.7 x 101	3.18	0.002
	log(Size Class)2	-1.8 x 101	5.0 x 100	-3.53	< 0.001
	Latitude	-4.9 x 104	2.3 x 104	-2.14	0.036
	Latitude2	6.4 x 102	3.0 x 102	2.14	0.036
	Depth	2.0 x 101	1.1 x 101	1.85	0.068
Diversity (Shannon H)	Intercept	1.4 x 103	7.9 x 102	1.72	0.090
	log(Size Class)	2.1 x 10-1	4.8 x 10-2	4.31	< 0.001
	log(Size Class)2	-3.7 x 10-2	8.9 x 10-3	-4.14	< 0.001
	Latitude	-7.1 x 101	4.1 x 101	-1.71	0.091
	Latitude2	9.2 x 10-1	5.4 x 10-1	1.71	0.091
	Depth	1.2 x 10-2	1.9 x 10-2	0.62	0.538
Evenness (Pielou J)	Intercept	-2.3 x 101	2.7 x 101	-0.87	0.386
	log(Size Class)	-4.0 x 10-3	1.6 x 10-3	-2.48	0.016
	log(Size Class)2	7.0 x 10-4	3.0 x 10-4	2.32	0.023
	Latitude	1.2 x 100	1.4 x 100	0.87	0.387
	Latitude2	-1.6 x 10-2	1.8 x 10-2	-0.87	0.388
	Depth	-3.0 x 10-4	6.4 x 10-4	-0.47	0.640

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013954 -0.005126 -0.002510  0.000872  0.105965 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.975e+01  2.757e+01  -0.716   0.4762  
log(Size_Class)      -3.283e-03  1.670e-03  -1.966   0.0533 .
I(log(Size_Class)^2)  5.746e-04  3.108e-04   1.849   0.0687 .
lat                   1.027e+00  1.436e+00   0.715   0.4771  
I(lat^2)             -1.334e-02  1.870e-02  -0.713   0.4780  
depth                -2.417e-04  6.615e-04  -0.365   0.7159  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0157 on 69 degrees of freedom
Multiple R-squared:  0.06698,   Adjusted R-squared:  -0.0006297 
F-statistic: 0.9907 on 5 and 69 DF,  p-value: 0.4299

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7 
 0.0116707844  0.0116707844  0.0135818312  0.0135818312  0.0133535000  0.0160100159  0.0160100159 
            8             9            10            11            12            13            14 
 0.0139468521  0.0098844429  0.0043185124  0.0043185124  0.0062295592  0.0062295592  0.0060012280 
           15            16            17            18            19            20            21 
 0.0086577439  0.0086577439  0.0065945801  0.0065945801  0.0011020710  0.0011020710  0.0030131178 
           22            23            24            25            26            27            28 
 0.0030131178  0.0027847867  0.0054413025  0.0054413025  0.0033781387 -0.0006842705 -0.0006842705 
           29            30            31            32            33            34            35 
 0.0002187110  0.0002187110  0.0021297578  0.0021297578  0.0019014267  0.0019014267  0.0045579425 
           36            37            38            39            40            41            42 
 0.0045579425  0.0024947787  0.0024947787  0.0009197555  0.0028308023  0.0028308023  0.0026024712 
           43            44            45            46            47            48            49 
 0.0026024712  0.0052589870  0.0052589870  0.0052589870  0.0031958232  0.0031958232 -0.0008665860 
           50            51            52            53            54            55            56 
-0.0008665860  0.0033429274  0.0033429274  0.0052539742  0.0052539742  0.0050256431  0.0050256431 
           57            58            59            60            61            62            63 
 0.0076821589  0.0076821589  0.0076821589  0.0056189951  0.0056189951  0.0015565859  0.0015565859 
           64            65            66            67            68            69            70 
 0.0066852276  0.0085962744  0.0085962744  0.0083679432  0.0083679432  0.0110244591  0.0110244591 
           71            72            73            74            75 
 0.0110244591  0.0089612953  0.0089612953  0.0048988861  0.0048988861 

$se.fit
 [1] 0.005493164 0.005493164 0.005343397 0.005343397 0.005314882 0.005654413 0.005654413 0.005971910
 [9] 0.006927905 0.003853380 0.003853380 0.003727927 0.003727927 0.003500958 0.004078885 0.004078885
[17] 0.004376629 0.004376629 0.003918866 0.003918866 0.003824869 0.003824869 0.003508152 0.004088390
[25] 0.004088390 0.004315988 0.005735008 0.005735008 0.003918423 0.003918423 0.003818133 0.003818133
[33] 0.003452327 0.003452327 0.003991513 0.003991513 0.004187817 0.004187817 0.003778022 0.003647954
[41] 0.003647954 0.003256201 0.003256201 0.003753828 0.003753828 0.003753828 0.003955533 0.003955533
[49] 0.005404797 0.005404797 0.003930715 0.003930715 0.003750184 0.003750184 0.003395068 0.003395068
[57] 0.003750491 0.003750491 0.003750491 0.003968722 0.003968722 0.005315321 0.005315321 0.004978438
[65] 0.004784446 0.004784446 0.004548023 0.004548023 0.004710012 0.004710012 0.004710012 0.004913433
[73] 0.004913433 0.005948707 0.005948707

$df
[1] 69

$residual.scale
[1] 0.015702

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

Metric	Term	Estimate	Standard Error	T Value	p
Richness (Breakaway – Betta)	Intercept	7.1 x 106	2.4 x 102	NA	< 0.001
	log(Size Class)	1.2 x 102	6.1 x 101	NA	0.058
	log(Size Class)2	-5.0 x 101	1.2 x 101	NA	< 0.001
	Latitude	-3.7 x 105	6.1 x 100	NA	< 0.001
	Latitude2	4.8 x 103	1.6 x 10-1	NA	< 0.001
	Depth	1.5 x 102	1.0 x 101	NA	< 0.001
Richness (Breakaway – LM)	Intercept	7.1 x 106	3.3 x 106	2.13	0.036
	log(Size Class)	2.4 x 102	2.0 x 102	1.21	0.231
	log(Size Class)2	-7.5 x 101	3.8 x 101	-2.00	0.050
	Latitude	-3.7 x 105	1.7 x 105	-2.13	0.037
	Latitude2	4.8 x 103	2.3 x 103	2.13	0.037
	Depth	1.5 x 102	8.0 x 101	1.89	0.064
Diversity (Shannon H)	Intercept	1.4 x 103	7.9 x 102	1.72	0.090
	log(Size Class)	2.1 x 10-1	4.8 x 10-2	4.31	< 0.001
	log(Size Class)2	-3.7 x 10-2	8.9 x 10-3	-4.14	< 0.001
	Latitude	-7.1 x 101	4.1 x 101	-1.71	0.091
	Latitude2	9.2 x 10-1	5.4 x 10-1	1.71	0.091
	Depth	1.2 x 10-2	1.9 x 10-2	0.62	0.538
Evenness (Pielou J)	Intercept	-2.0 x 101	2.8 x 101	-0.72	0.476
	log(Size Class)	-3.3 x 10-3	1.7 x 10-3	-1.97	0.053
	log(Size Class)2	5.7 x 10-4	3.1 x 10-4	1.85	0.069
	Latitude	1.0 x 100	1.4 x 100	0.71	0.477
	Latitude2	-1.3 x 10-2	1.9 x 10-2	-0.71	0.478
	Depth	-2.4 x 10-4	6.6 x 10-4	-0.37	0.716

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)